Helmholtz Problem

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *

keys: Acoustic Scattering, Brakhage-Werner formulation, Helmholtz equation, Dirichlet data

Helmholtz Problem#

We consider a Helmholtz problem with a plane wave as initial condition, i.e., \(m = e^{i\,\kappa \,x} \) and

\[\begin{split} \left\{ \begin{array}{rcl l} \Delta u + \kappa^2 \, u &=&0, \quad &\Omega \subset \mathbb R^3\,,\\ \gamma_0\, u &=& m, \quad &\Gamma = \partial \Omega\, .\end{array} \right. \end{split}\]

Combined field integral equations combine single and double layer integral operators, one simple option is the Brakhage-Werner formulation.

The solution is represented as

\[ u = i \,\kappa \,\mathrm{SL}(j) - \mathrm{DL}(j)\,, \]

where \(j\) solve the boundary integral equation $\( \left( \frac{1}{2} + \mathrm K + i \, \kappa \,\mathrm V \right) \, \mathrm j = \mathrm u_{in} \qquad \text{on} \, \Gamma = \partial \Omega \)$

The single layer and double layer operator are defined using the fundamental solution of the Helmholtz operator,

\[ G(x,y) = \frac{e^{i\kappa \|x-y\|}}{4\pi\,\|x-y\|}. \]

The layer potentials \( \mathrm SL\) and \(\mathrm DL\) have the same form as in the Laplace case, except that the Laplace kernel is replaced by the Helmholtz kernel \(G\).

kappa=10
order=4
screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()

sp = Fuse(Sphere( (0,0,0), pi).faces)
screen.faces.name="screen"
sp.faces.name="sphere"
shape = Compound([screen,sp])

mesh = shape.GenerateMesh(maxh=5/kappa).Curve(order)
Draw (mesh);
fes_sphere = Compress(SurfaceL2(mesh, 
                                order=order, 
                                complex=True, 
                                definedon=mesh.Boundaries("sphere")))
u,v = fes_sphere.TnT()
fes_screen = Compress(SurfaceL2(mesh, order=order, 
                                dual_mapping=True, 
                                complex=True, 
                                definedon=mesh.Boundaries("screen")))
print ("ndof_sphere = ", fes_sphere.ndof, "ndof_screen =", fes_screen.ndof)
ndof_sphere =  17010 ndof_screen = 31110
with TaskManager():
    C = HelmholtzCF(u*ds("sphere"), kappa)*v*ds
    u,v  = fes_sphere.TnT()
    Id = BilinearForm(u*v*ds).Assemble()
lhs = 0.5 * Id.mat + C.mat
m = exp(1j * kappa * x) 
rhs = LinearForm(m*v*ds).Assemble()
gfu = GridFunction(fes_sphere)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    gfu.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=40, tol=1e-8)
GMRes iteration 1, residual = 60.52440579389444     
GMRes iteration 2, residual = 24.001710415059282     
GMRes iteration 3, residual = 10.81618663970256     
GMRes iteration 4, residual = 5.17850130118241     
GMRes iteration 5, residual = 2.549268448878951     
GMRes iteration 6, residual = 1.288805495688266     
GMRes iteration 7, residual = 0.6549832602491101     
GMRes iteration 8, residual = 0.3554164558330166     
GMRes iteration 9, residual = 0.19990866606606392     
GMRes iteration 10, residual = 0.10836869993844735     
GMRes iteration 11, residual = 0.059314453140396746     
GMRes iteration 12, residual = 0.032946081231833756     
GMRes iteration 13, residual = 0.013819139108246862     
GMRes iteration 14, residual = 0.007360642470274827     
GMRes iteration 15, residual = 0.0033138105281745843     
GMRes iteration 16, residual = 0.0019472886266694725     
GMRes iteration 17, residual = 0.001152960819267155     
GMRes iteration 18, residual = 0.0006865451125338632     
GMRes iteration 19, residual = 0.00040647668374400083     
GMRes iteration 20, residual = 0.00023525302450811517     
GMRes iteration 21, residual = 0.00012568335128557858     
GMRes iteration 22, residual = 7.27284483870596e-05     
GMRes iteration 23, residual = 4.282868528421723e-05     
GMRes iteration 24, residual = 2.238947095910078e-05     
GMRes iteration 25, residual = 1.1323654081077384e-05     
GMRes iteration 26, residual = 5.6817723197032995e-06     
GMRes iteration 27, residual = 2.815481215969824e-06     
GMRes iteration 28, residual = 1.3918984915180687e-06     
GMRes iteration 29, residual = 7.931632776571698e-07     
GMRes iteration 30, residual = 4.523061248529819e-07     
GMRes iteration 31, residual = 2.723734968721705e-07     
GMRes iteration 32, residual = 1.6096117668515718e-07     
GMRes iteration 33, residual = 9.476709863740349e-08     
GMRes iteration 34, residual = 5.214502615501129e-08     
GMRes iteration 35, residual = 2.9709338426800344e-08     
GMRes iteration 36, residual = 1.733994073439006e-08     
GMRes iteration 37, residual = 9.167410459121383e-09     
Draw (gfu, order=5, min=-1, max=1);

Postprocessing on screen#

uscat = GridFunction(fes_screen)
with TaskManager():
    uscat.Set(1j*kappa*HelmholtzSL(u*ds("sphere"),kappa)(gfu) - HelmholtzDL(u*ds("sphere"), kappa)(gfu), \
              definedon=mesh.Boundaries("screen"))
print ("Scattered field")
Draw (uscat, mesh, min=-1,max=1, animate_complex=True, order=4);
Scattered field
uin = mesh.BoundaryCF( {"screen": m }, default=0)
print ("Total field")
Draw (uin+uscat, mesh, min=-1,max=1, animate_complex=True, order=4);
Total field